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Abstract 

Surnames and nonrecombining alleles are inherited from a single parent 
in a highly similar way. A simple birth-death model with mutations can 
accurately describe this process. Exponentially growing and constant pop- 
ulations are investigated, and we study how different compositions of the 
founder populations can be observed in present-day diversity distributions. 
We analyse different quantities in the statistically stationary state, both 
through analytic and numerical methods. Our results compare favourably 
to field data for family sizes in several countries. We discuss the relation- 
ship between the distribution of surnames and the genetic diversity of a 
population. 

1 Introduction 

Biological and cultural features of human populations have been traditionally 
studied by separate disciplines, but the parallelisms between biological and cul- 
tural evolution have been put forward by a number of researchers. Already Dar- 
win (1871) pointed out that "the formation of different languages and of distinct 
species and the proofs that both have been developed through a gradual process are 
curiously parallel. " 

Cultural traits are transmitted from ancestors to their descendence, in a pro- 
cess analogous to inheritance, and are subject to changes, similar to mutations, by 
interaction between individuals -such as teaching and imitation. Moreover, they 
usually fulfill a practical purpose, which amounts to being subject to selection. In 
fact, they enhance the relationships within human groups, defining social entities 



1 



comparable to certain biological species and populations. The quantification of 
cultural traits has been attempted only recently. For example, Cavalli-Sforza et 
al. (1982) applied concepts from the quantitative theory of biological evolution 
to construct a theory of cultural evolution. They analyzed forty traits, ranging 
from political preferences to superstitions. Many of these traits are subject to 
high mutability, since they are influenced by fashion, individual acquaintances, 
and personal experience, and one does not expect their quantitative properties 
to be directly comparable to those of any biological feature. Other traits, on 
the other hand, are better preserved. Among them we find languages and sur- 
names. Language is essential to integrate the individual to society; surnames are 
historical -though recent- signs of identity in social groups. Quite early, Galton 
and Watson (1874) dealt with the problem of the extinction of surnames. This 
problem is equivalent to that of the extinction of a mutant allele in a population, 
although this relation was only established half a century later (Fisher, 1922; 
Haldane, 1927), when the first quantitative approaches to biological evolutionary 
processes took place. 

Comparative methods analogous to biological taxonomy are used to determine 
the degree of similarity between languages. This method returns the genetic clas- 
sification of linguistic diversity (sec for example Grecnberg, 1992; Ruhlcn, 1992). 
Recently, a quantitative study of the taxonomy of languages has been carried out 
(Zanette, 2001), showing that the distribution of the number of subtaxa within 
a taxa displays scaling properties, quantitatively similar to those disclosed in 
biological taxonomy by Burlando (1990, 1993). That is, if n is the number of 
subtaxa belonging to a given taxa -say, the number of languages in the Indo- 
European family, or the number of species in the genus Canis- the fraction p{n) 
of taxa that have precisely n subtaxa scales as p{n) ~ n~^. The exponent is 
found in the interval 1 < /3 < 2 in both cases. This gives quantitative support 
to Darwin's observation on the "equivalence" between the mechanisms behind 
biological and linguistic evolution. A complementary comparison is that of lin- 
guistic abundance, measured as the number of individuals speaking a language, 
with the number of individuals of a biological species. Again, the frequency as a 
function of the number of individuals has scaling properties both for languages 
(Gomes et al., 1999) and for species (Pielou, 1969), see Fig. 1. 

Surnames arc cultural traits (Cavalli-Sforza & Feldman, 1981) whose trans- 
mission bears strong similarity with that of some biological features. They are 
inherited from one of the parents, usually the father, much in the same way as 
the Y-chromosome or the mitochondrial DNA. The extinction of a surname and 
the persistence of a non-recombining neutral allele are equivalent problems. This 
is not only a mathematical fact, but has also practical implications. Indeed, to 
assess the multiple or single origin of a surname one can turn to genetic measures. 
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Figure 1: Scaling behaviour of the fraction p{n) of species represented by n 
individuals (triangles), of surnames borne by n persons (squares), and of lan- 
guages with n speakers (circles). Data for the distribution of species abun- 
dance are from Poore (1968), corresponding to trees in a Malaysian rain- 
forest [see also Sole and Alonso (1998)]. Data for surnames (beginning by 
A) are from the 1996 Berlin telephone book. Data for languages are from 
|http : //www, sil . org/ethnologue/pref ace .html| . As a guide to the eye we 
draw lines with slope —2 (solid line) and —3/2 (dotted lines). Data sets have 
been mutually shifted for better visualization. 

since males sharing the same surname might also share the same haplotype in the 
nonrecombining segment of the Y chromosome (Sykes & Irven, 2000). In a very 
large population, the statistical properties of the surname distribution can be 
strongly correlated with genetic diversity (Barrai et al, 1996), and may even be 
used to understand the genetic structure of a population (Yasuda et al., 1974). 
Recent reports on actual populations (Miyazima et al, 2000; Zanette &: Man- 
rubia, 2001) show that the distribution of surnames follows the same statistical 
law observed for languages and biological species. Namely, if n is the number of 
individuals bearing a given surname, the fraction p{n) of surnames decreases with 
n as p(n) ~ n~'^ (see Fig. 1). Here, however, the exponent is always 7 ~ 2. This 
paper focuses on a model aimed at predicting this kind of regularities, observed 
in disparate human populations. 



3 



2 Theoretical approaches to surname evolution 



At the mathematical level, several models have been proposed and analysed in 
order to identify quantitative properties of surnames evolution -or, equivalently, 
those of nonrecombining neutral alleles. The main points addressed by these 
studies are (i) the probability of fixation of a given surname/allele in a closed 
population (Galton & Watson, 1874; Fisher, 1922; Haldane, 1927; Moran, 1962; 
Lange, 1981; Rannala 1997; Hull 1998), and (ii) the distribution of the number 
of individuals bearing the same surname/allele (Kimura & Crow, 1964; Karlin 
& McGregor, 1967; Fox & Lasker, 1983; Panaretos, 1989; Gale 1990; Consul, 
1991; Islam, 1995; Zanette & Manrubia, 2001). Indeed, these two questions cover 
complementary aspects of the same problem. In (i), one deals with a closed 
population (no immigrants enter the system), and implicitly assumes that the 
mutation rate is low enough, such that fixation can indeed occur before a mutant 
form appears. Suppose that there are N individuals in the population. We know 
from coalescence theory that the time Ug (in units of the number of generations) 
required for a neutral allele to be fixed is of order Ug ^ N. Now suppose that 
the mutation rate per generation and per individual is r. Then M = rUgN is the 
average number of mutants after rig generations have elapsed. Only if M <^ 1, 
that is if r ^ N"^, will the fixation of the allele be possible. As soon as this 
inequality is violated, a new situation arises, in which both neutral drift and 
mutation play relevant roles. In this case, a broad distribution of surnames or 
of genetic diversity arc expected. Actually, the value of r will be usually fixed 
by the nature of the problem, while the size N of the considered population can 
increase enough such that > 1/r. Therefore, the statistical behaviour crosses 
over to the second regime, where the appearance of mutants cannot be discarded. 
We are then in the assumptions of the models of class (ii). 

All the models in class (ii) represent systems where statistical equilibrium 
is reached. While the neutral drift drives the less frequent surnames out of 
the system, mutations generate new surnames. For large times, the number 
of different surnames in the system is much larger than unity. Yasuda and co- 
workers (Yasuda et al, 1974) used a stochastic model for a population with a fixed 
number of individuals. At each evolution step, a new individual with the surname 
of his father is added, but he acquires a new surname with a given probability. 
In any case, a randomly chosen individual is eliminated in order to keep the 
population constant. Their analytical results compared successfully to field data 
despite the restriction of constant population, which forces the elimination of 
individuals. In this model, the size of the progeny is Poisson distributed. 

Other models start with a single individual in the population and new individ- 
uals carrying the same surname or a new one are sequentially added (Panaretos, 
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1989; Consul, 1991; Islam, 1995). These fall in the category of branching pro- 
cesses (Harris, 1963) with an increasing number of individuals in the population. 
Panaretos (1989) rephrased a model introduced by Simon (1955) in the context 
of linguistics, as follows. At each step, an individual is chosen to be the father 
of a newborn and his surname is transmitted with probability 1 — a. With the 
complementary probability a, the newborn is assigned a surname not present in 
the population. Recently, we have modified Simon's model by (i) introducing 
an additional parameter // which represents the death rate of the individuals in 
the population, and (ii) allowing for arbitrary initial conditions. Preliminary re- 
sults have been successfully compared with large data sets (Zanette & Manrubia, 
2001). This is the starting point of the present work, where we present new 
analytical and numerical results for this birth-death model. 

We describe the model in the next section, where our analytical results are 
derived. Questions like the role played by the composition of the founder popula- 
tion and by death events in an exponentially growing population are addressed. 
We also analyse a system with constant population, which turns out to quan- 
titatively differ from the previous case. In Section ^ we numerically check our 
analytical results and run computer simulations of the model in situations not 
amenable to an analytical description. In Section |5| our results are tested against 
field data. We finish with an overall discussion in the last section. 



3 Birth-death model for surname evolution 

Our model population evolves in discrete steps, each step corresponding to the 
birth of a new individual. At each time step, moreover, an individual is chosen 
at random from the whole population and becomes removed with probability /i, 
representing a death event. The total population at step s, P{s), is therefore a 
stochastic process governed by the evolution equation 

Pis + l) = Pis) + l-wis), (1) 

where w{s) is the dichotomic stochastic process 

, . J 1 with probability fi, , . 

jo with probability 1 — /i. 

We show in Appendix ^ that, averaging over realizations of the stochastic process 
w{s), the average total population P grows exponentially in time: 

P{t) = Poexp[i^il-fi)t]. (3) 

Here, u is the birth rate per individual and unit time, and the product i^fi turns 
out to be the corresponding death rate; Pq is the initial population. 
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We think of the population as divided into groups -the famihes- within which 
ah individuals bear the same surname. At each birth event, the newborn is as- 
signed a new surname, not previously present in the population, with probability 
a. This probability can be seen as a mutation rate for surnames, but could also 
be interpreted as a measure of migration towards the population of individuals 
with new surnames. With the complementary probability, 1 — a, a preexistent 
individual is chosen at random from the population and becomes the newborn's 
father, i.e. his surname is given to the newborn. Surnames are therefore assigned 
with probabilities proportional to the size of the corresponding families, allow- 
ing however for fluctuations inherent to the random distribution of births among 
families. Consequently, the distribution of surnames is driven by a stochastic 
multiplicative process (Van Kampen, 1981), modulated in turn by the total pop- 
ulation growth. This process is analogous to the mechanism proposed by Simon 
(1955) to model the frequency distribution of words and city sizes, described by 
Zipf's law (Zipf, 1949), among other instances. Our model, in fact, reduces to 
Simon's model if mortality is neglected, i.e. for /i = 0. Note that allowing for 
death events adds a relevant process in the case of surnames, namely, the pos- 
sibility that a surname disappears if it is borne by a single individual and the 
individual dies. 

Since neither the probability of becoming father of a newborn nor the death 
probability depend on the individual's age, the present model population can be 
thought of as ageless. As a consequence, if and fi are constant, the probability 
that an individual dies at an age between T and T -|- dT is 



which implies a life expectancy T = Moreover, the probability that an 

individual has exactly m children with the same family name during its whole 
life equals 



giving a fertility m = 1/(1 — a)/x. The exponential distribution in (|5|) is in 
reasonable agreement with actual data collected over relatively short periods, for 
instance, for the United States (Hull, 1998), but contrasts with the Poissonian 
distribution found in data integrated over several centuries, as in the case of 
England (Dewdney, 1986). This discrepancy can presumably be ascribed to long- 
range variations of real birth and death rates (see Fig. 2). 

We are assuming that an individual's surname is inherited from the father. 
Consequently, the model -as presented above- describes the evolution of the male 
population only. The same evolution rules apply however if it is assumed that 
surnames are transmitted with the same probability by either parent. In this 
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number of children 

Figure 2: Probability distribution of the number of children per male in two 
different populations. Open dots correspond to an average during the period 
1350-1986 in England (Dewdney, 1986). The solid line joining the data points is 
actually a fitting with a Poisson distribution, P{n) = exp(— A)A"/n! with average 
(n) = A = 1.15. Solid dots arc data from the 1920 American census [from Hull 
(1998); originally compiled by Lotka (1931)]. The solid line fitting the first part 
of the data is an exponential distribution, P{n) oc exp(— 0.3n). 

case, there is no sex distinction and the model encompasses the whole popula- 
tion. The real situation is in fact intermediate between these two limiting cases: 
whereas some societies strongly favor inheritance of the surname cither from the 
father or from the mother, other groups allow for a choice between the two possi- 
bilities. In some developed western countries, where the surname is preferentially 
transmitted by the father, opting for the mother's surname has been proposed 
as a method to avoid the persistent extinction of surnames, due to the extremely 
slow population growth (Legay & Vernay, 1999). 

3.1 Distribution of families by size. The case of = 

As stated above, our model reduces to Simon's model (Simon, 1955) for /j, = 0. 
In this case, u)(s) = for all s and the total population evolves deterministically, 
P(s) = Pq + s, since exactly one individual is added to the population at each 
step. The number of families with exactly i individuals at step s, ni{s), grows to 
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ni{s) + 1 when an individual is added to a family of size i — 1. This occurs with 
probability (1 — a){i — 1)/P(s). On the average, thus, the number of families 
withi individuals varies according to 

_ _ I — a _ _ 

ni{s + 1) = ni{s) + -^j-j [{i - l)ni_i(s) - mj(s)] , (6) 

for i > 1. To the families with only one individual, on the other hand, the 
positive contribution comes from the creation of new surnames with probability 
a. Therefore, 

Note that, since for fi = the stochastic process w{s) becomes trivial, overlines 
indicate here average over realizations of the stochastic process by which each 
newborn's surname is chosen. 

The system of equations (^) and (0) can be completely solved for arbitrary 
initial conditions. In fact, (^) is an autonomous equation for ni(s), whose solution 
reads 

ni{s) = Po + s)+( m - Po p p ^\wp 7—4 , 8 

2 — a V 2 — a / L (Pq + s)l (Pq — 1 + a) 

where T{z) is the gamma- function (Abramowitz & Stegun, 1970). Then, eqn (|^ 
can be used to recursively obtain nj(s) for i > 1. For long times, s — > oo, rli(s) is 
essentially given by the first term in the right-hand side of eqn @, which grows 
linearly with s, as the total population. The second term, which contains all 
the information about the initial condition, can be seen to decrease as s~^^~°'\ 
In this limit, the recursion from eqn (^) can be immediately solved. Thus, for 
s — > oo, we find 

Ms) = V> ( (^0 + s). (9) 



a 



The asymptotic number of families of a given size i turns out to be proportional to 
the total population. For fixed s and sufficiently large values of i, nj(s) decreases 
as a power law, flj oc with 

z = l + -^. (10) 

1 — a 

In the limit of small mutation rate, a ~ 0, we have z ~ 2. As advanced in the In- 
troduction, this is the value observed in actual surname distributions. This same 
result was obtained by Simon (1955) for a special initial condition. Elsewhere 
(Zanette &: Manrubia, 2001) we have already discussed the fact that, though the 
effects of the initial condition in Simon's model fade out for long times, transients 
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can strongly depend on the initial distribution nj(0) (see also section |^. In our 
case, this aspect could be relevant in populations where modern surnames have 
appeared recently, like the Japanese (Miyazima et al, 2000). We shall return to 
this specific point later. 

Equations @ and (0) imply that, in average, the total number of surnames 
in the population grows linearly: 

oo 

N{s) = J2ni{s) = No + as, (11) 

1=1 

with A^o the initial number of different surnames. As a function of time, the 
number of surnames increases exponentially: 

N{t) = No + aPo [exp(zyt) - 1] . (12) 

We have already pointed out that, in contrast, the number of surnames in some 
real populations at present times is decreasing (Cavalli-Sforza & Feldman, 1981; 
Legay &: Vernay, 1999). Under suitable conditions, adding mortality allows to 
reproduce this particular behaviour. 



3.2 Effect of death events 

Under the action of mortality, the growth of the total population P{s) fluctuates 
stochastically, according to eqn (||) , depending on the occurrence of death events 
at each evolution step. The evolution of nj(s) can be implemented in two sub- 
steps, as follows. First, eqns and (|^ are used to calculate the intermediate 
values 

_ _ I — a _ _ 

n'i{s) = riiis) + [{i - l)ni_i(s) - ini{s)] , (13) 

and 

n' {s) = ni{s)+a- -^r-rUi (s) , (14) 

and the population is updated to P'{s) = P{s) + 1. Second, the evolution due to 
mortality is applied. In terms of the random process w{s) of eqn (|2|), we have 



[{i + l)n[^^{s)-in[{s)], (15) 

and P{s + 1) = P'{s) - w{s) [cf. eqn (jl])]. 

Naive replacement of the stochastic process w{s) by its average value, w{s) 
H, in order to obtain a deterministic approximation to the problem, would lead 
to an equation for which positive solutions cannot be insured. In fact, it is 
possible to give initial conditions for such deterministic equation which lead to 



nAs + 1) 



■w[s) 



P'(s) 
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negative values of nj(s) at sufficiently large s and i. A satisfactory deterministic 
approximation can however be proposed by assuming that the solution ni(s) 
varies slowly both in s and i. This assumption is generically verified for large 
populations with smooth initial conditions. In this situation, eqns (^) and ( |T5[) 
admit a continuous approximation in terms of real variables z and y, which replace 
the integer variables s and i, respectively. Meanwhile, ni(s) is replaced by a 
continuous function n{y,z). 

The continuous approximation can be analyzed at different truncation orders, 
as discussed in Appendix p|. At the first order, the approximate solution to eqns 



(y) and (|15D reads 



for y < yriz), and 



a — iJ. 



n{y, z) = yj,^n{y/yT{z), 0) (17) 



1 _ ,y \ {1-Q!-/')/(1-M) 

yT(^)=(l + — . (18) 



for y > yriz), with 



In this first-order continuous approximation, thus, the average number of fam- 
ilies n{y, z) exhibits two separated regimes. For large values of the size variable, 
y > yriz), the distribution is essentially determined by the initial condition. At 
the first evolution steps, where z — > and yr ^ I, this regime covers practically 
all the domain of variable y. As time elapses and yx grows, however, this regime 
recedes and is replaced for y < yx hy & power-law distribution, n a y~^, with 

C = l+ , ^"^ • (19) 

1 — a — fj. 

Note that, as the exponent z in eqn ([lO|), this new exponent coincides with the 
observed value, ~ 2, in the relevant limit a « 0. This result is independent 
of the death probability ii. For sufficiently long times, the power-law regime will 
be observed for all the relevant family sizes. All the contribution of the initial 
condition will become restricted to the domain of the largest families. 

The boundary yriz) between the two regimes, eqn (|T^), grows as a power of 
the ratio between the average current population Pq + (1 — fj.)z and the initial 
population Pq. For a 0, the exponent is practically equal to unity. Conse- 
quently, for the boundary to reach a given value yo, in such a way that all the 
families with sizes below yo are in the power-law regime, the total population 
must grow by a factor practically equal to yo- Taking into account eqn (|3|), the 
position of the boundary as a function of real time reads 

yT(t) =exp[i/(l -a-/x)t]. (20) 
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The transient to needed for the power-law regime to develop up to a given size 
j/o is therefore logarithmic, oc Inyo- 

The average total number of surnames in the continuous approximation is 
calculated by analogy with eqn (|Tll): 

/oo 
n{y,z)dy. (21) 

Using the above first-order solution for n{y,z), we find N{z) = Nq + az or, as a 
function of time, 

N{t) = No + -^{exp[z.(l - f,)t] - 1}, (22) 
i fi 

as shown in Appendix ^ [cf. eqn ([l^)]. In consequence, within this approxima- 
tion, the surname diversity always grows exponentially. We point out, however, 
that this conclusion is valid for sufficiently smooth distributions and for fi < 1 — a, 
two necessary conditions for the continuous approximation to apply to our prob- 
lem. It could therefore be argued that having a decreasing number of surnames, 
as found in some modern developed societies (Cavalli-Sforza &: Feldman, 1981; 
Legay & Vernay, 1999), requires violation of the above conditions. For instance, 
an initial condition with sharp variations would violate the smoothness condition 
during the first evolution stages. A death probability ^ ~ 1 would also threaten 
the validity of the continuous approximation. This is precisely the case of mod- 
ern developed societies, where the population growth rate is practically vanishing 
and, as a matter of fact, reaches occasionally negative values -a situation not de- 
scribed by the present model. In Section we treat the special case = 1 and 
show that, in this limit, the total number of surnames can in fact decrease. 

Though the first-order continuous approximation gives a rather rough de- 
scription of the solution to our model, as a piecewise function with a moving 
discontinuity, it provides a quite clear qualitative picture of how the solution be- 
haves. The growth of the power-law regime as the initial condition recedes is in 
good qualitative agreement with the evolution observed in numerical realizations 
of the model. A better analytical approximation is obtained from the second- 
order truncation. Second-order derivatives would in fact introduce diffusive-like 
effects in the variable y, with the consequent smoothing of the discontinuities of 
the first-order approximation. The second-order equation, however, cannot be 
analytically solved for arbitrary initial conditions. Nevertheless, it is possible to 
give an asymptotic approximation for long times, as 

aP(z) ( 1 — a — 1 — a — yu\ 

n[y,z) = - ^ 2- y-if/ C-1,0,2- -^y , (23) 

\ — a — — a^[i) \ 1 — a + / 

where U{a,b,x) is the Kummer's function (see Appendix P). With respect to 
the first-order approximation, eqn ([l6|), this solution predicts a lower value of 
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n{y, z) for small y. For larger family sizes, however, it behaves as a power law, 
n{y,z) oc y~^ with exactly the same exponent as in eqn (|T9|). The asymptotic 
behaviour for large y in the power-law regime is therefore not modified. It must 



be pointed out that, for small values of a, the profile of n{y, z) given by eqn (23) 
results to be quite independent of the death probability over a considerable range 
of values of /i. Only for /x ~ 1 — a ~ 1, where the exponent strongly depends on 
/i, does n{y, z) change sensibly as fi is varied. We take advantage of this feature 
in Section ^ where the second-order approximation is compared with actual data 
of surname distributions. 

3.3 The case of constant population, /i = 1 

As argued above, the limit /i = 1 is relevant to the discussion of the evolution 
of surname distributions in some modern developed societies. In fact, this limit 
corresponds to the case where the population growth rate vanishes, a situation 
which is closely approached, for instance, in many European countries. In this 



case, eqn (15) is again deterministic, with w{s) = 1 and P'{s) = Pq + 1. The 
total population at any time is P{s) = Pq. 

The main difference between this case and that of /i < 1 is that, now, the 
distribution nj(s) becomes independent of time for asymptotically large times. 
This feature is in agreement with the fact that the asymptotic distribution for 
/X < 1 is proportional to the total population P{s). The asymptotic surname 
distribution 7t?° is given, for /i = 0, by the recurrence equations 



_oo ^ [(2 - a)Po - 2i{l - a)]inr - (1 - a)(Po + 1 - - l)nti ..4^ 
"^+^ (i + l)[Po-(l-«)(i + l)] ^ ' 

for i > 1, and 



-oo _ [(2-a)Po-2(l-a)]nr-aPo^ 

~ 2[Po-2(l-a)(z + l)] ■ ^^^^ 

Note that, since the total population is always Pq, we have = for i > Pqi 
as no family can be larger than Pq. 

The solution reveals two well-defined regimes, depending on how the product 
aPo compares with unity. For aPo > 1, the asymptotic distribution behaves as 

nr^^il-ay-' (26) 

for a vast range of family sizes. Departures from this behaviour are found very 
close to i = Pq only. For i > Pq, in fact, the distribution must vanish. We stress 
the remarkable difference between the exponential stationary distribution ( p6|) 
and the long-time power-law solution obtained for // < 1. In this regime, the 
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stationary total number of surnames is 

N^^^\l^a\. (27) 
I — a 

For aPo < 1, on the other hand, the distribution behaves as a power law, 
~ over practically the whole range of family sizes. Note that the exponent 
of this power law is in agreement with eqn (|T^ ) in the limit /i ^ 1. For i ~ Pq; 
however, the distribution deviates from the power law and exhibits a sharp peak. 
In the limit a — > 0, the distribution becomes an isolated peak at i = Pq, namely, 
n?° = for i < Pq and nf^ = 1 for i = Pq. Therefore, the total number 
of surnames is N°° = 1. This special solution describes the well-known case 
of a closed population with no surname mutations where, by random drift, a 
single surname survives for asymptotically large times (Cavalli-Sforza & Feldman, 
1981). Notice also that this limit corresponds to the birth-death model introduced 
by Moran (1962) to study probabilities of fixation of alleles when generations 
overlap. 

We conclude that, for a given population Pq, the asymptotic number of sur- 
names can be very small if the mutation probability a is, in turn, small enough. 
Consequently, in a steady population with many surnames at the initial stage 
and with a sufficiently low mutation probability, the number of surnames will 
decrease towards the stationary value as time elapses. Numerical simulations, 
discussed in detail in the next section, show that for /i just below unity and small 
a, an initial transient where the number of surnames decreases temporarily can 
be observed. Since fi < 1, however, the population grows steadily and, as a 
consequence of mutations, will also increase in the long run. The situation 
in modern developed societies is that the population growth rate has been con- 
stantly decreasing, to reach values around zero at present times. Starting from a 
state with a wealth of surnames -due to the combined effect, a few centuries ago, 
of a high population growth and the frequent appearance of new surnames- these 
societies have now reached a regime of almost stationary population where the 
number of surnames decreases. This situation would be reverted if the growth 
rate could be maintained above zero during substantial periods. Presently, the 
only mechanism acting in this direction seems to be immigration, which is in 
addition an effective source of new surnames. 

4 Numerical results 

In this section we present results of numerical realizations of our model, in order 



to compare with the analytical approximations presented in Section 3.2, and 



to illustrate the behaviour within the regimes where such approximations do 
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not hold. Emphasis is put on the role of the initial conditions, the duration 
of transients, and the evolution of the number of different surnames. Also, we 
discuss actual surname distributions of modern populations in the light of our 
results. 



4.1 Role of initial conditions 

First, it is worthwhile to illustrate how the shape of the distribution nj(s) depends 
on the initial condition nj(0) [see also Zanette &: Manrubia (2001)]. We focus 
the attention on initial distributions of the form nj(0) = Nq for i = iq, and 
nj(0) = for i ^ zq, in which there are Nq families of equal size formed by 
iQ individuals each. Consequently, the initial population is Pq = ioNq. In the 
following, we denote such an initial condition by the pair (ig, A^o)- Figure 3 shows 
four normalized distributions of family sizes, 

obtained both from the numerical realization of the stochastic birth-death model 
and from the iterative solution of the deterministic eqns (^) and (^), for /i = 0. 
We find a very good agreement between both methods and, at the same time, 
clearly realize the relevant role of the initial condition in determining the profile 
of the distribution for large family sizes. These same features are found for other 
values of the death probability. 

The solution of the first-order approximation to our model, eqns ( [T^ ) to (|l8|), 
predicts that in the zone of small family sizes -i.e. for y < yx m the continu- 
ous variables- the only dependence on the initial condition appears through the 
quantity Pq. This means that the distribution of family sizes is not sensitive to 
a variation of and Nq as far as their product is kept constant. Note that the 
three cases (4, 5), (1, 20), and (20, 1) of Fig. 3 share the same value of Pq = 20. 
For the parameters of the figure, the crossover between the regions of small and 
large family sizes, given for the continuous variables by eqn (p^), should occur 
at It ~ 470. Though two of the distributions indeed have the same profile -with 
the expected power-law decay- up to that value, the distribution corresponding 
to the initial condition (1, 20) deviates considerably below it- This deviation can 
be attributed to the fact that the initial condition (1,20) corresponds to a quite 
singular distribution, with a high peak at i = 1. A continuous approximation 
for such a distribution is arguably expected to give a poor description of the real 
situation. 
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Figure 3: Effect of the initial condition on the tail of the surname distribution. 
The normalized distribution pi{s) is shown for s = 10^, // = 0, a = 10~^, and 
four different initial conditions, as specified in the legend. Symbols stand for 
numerical simulations of the birth-death stochastic model, averaged over 2000 
independent realizations. Curves correspond to the numerical solution of the 
deterministic equations (^) and (^). 



4.2 Computational measurement of transients 



Equation ( |18| ) defines, within the first-order approximation, the boundary that 
separates the asymptotic regime and the zone dominated by the initial condition. 
Alternatively, for a fixed family size, it can be used to determine the transient 
before the asymptotic distribution builds up at that size. In order to test this 
aspect of the first-order approximation, we have devised an independent com- 
putational method to evaluate the point at which the cut-off between the two 
regimes actually takes place. According to eqns ( |l^ ) and (|28|), the asymptotic 
normalized distribution of family sizes is p-^ = pi{s — > oo) — In 

numerical realizations, the distribution is expected to adopt values very close to 
p'^° for large i and large enough s. In the simulations, we fix a certain family 
size It and stop the calculation at the step st when the measured value of pi^, 
satisfies, for the first time, \pij. —p^\ < A. The results presented in the following 



correspond to the choice = 10^ and A = 10^^. Since, as shown in Section 4T 
the agreement between the numerical simulations of the stochastic birth-death 
model and the iterative solution of the corresponding average equations is very 
good, we use this second description to measure the transient st- 

We focus the analysis on the dependence with the initial conditions, and 
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keep the values of a = 10 ^ and /x = fixed. As far as a ~ 0, the results 
are qualitatively the same for other values of /i. We consider initial conditions 



of the form {io,NQ), as defined in Section 4.1. Figure 4 shows the transient 
St measured from the numerical solution to eqns and (0) according to the 
criterion introduced above, as a function of iq and Nq. On the one hand, for 
small A'o and, especially, for small io, st varies quite irregularly. On the other 
hand, we see that for larger values of io and A^'o, say io > 10 and A'^o > 5, st 
exhibits a well-defined linear dependence with both parameters. In this regime, 
the transient is well approximated by 

ST = aioNo. (29) 



Linear fitting of st as a function of io and Nq for the data shown in Fig. 4 
yields a = 237 it 4. Note that eqn ( p9| ) implies that, for large zq and Nq the 
transient length is proportional to the initial population Pq = ioNq. This is to 
be compared with the analogous result within the first-order approximation, eqn 
(p^), which implies 

(30) 

The first-order approximation predicts then a proportionality between st and Pq 



ST 



l/rp 1. 
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-independently of the values of a and //- in full agreement with the numerical re- 
sults. We point out, however, that the corresponding proportionality coefficients 
cannot be directly compared. In fact, the coefficient a in eqn ( p9D is evaluated 
following the computational definition of the transient, and thus depends on the 
parameter A, which has no correlate in the analytical approach. Note, neverthe- 
less, that for a and fi = eqn ( |30[ ) predicts a coefficient approximately equal 
to It = 100, which is of the same order as the numerical result. 



4.3 Evolution of surname diversity 

As discussed in Section |3|, the evolution of the average number N{s) of different 
surnames is driven in our model by two competing mechanisms. The diversity 
increases due to the appearance of new surnames at rate a, and surnames borne 
by single individuals disappear when the individual in question dies. Within the 
first-order continuous approximation to our model, for /i < 1, the average number 
of different surnames increases linearly with s as N(s) = Nq + as (see Section 
|3.2| ). In this approximation, the death probability /j plays a role in the variation 
of diversity as a function of time, eqn (p^. On the one hand, increase of the 
surname diversity in a steadily growing population is generally expected when 
new surnames are created at a constant rate. On the other, it is possible to 
conceive special situations where the number of different surnames should tem- 
porarily decrease, violating the first-order approximation. Imagine, for instance, 
that the initial population consists of an ensemble of families with only one in- 
dividual each. For moderate values of /x and small a, the initial stage would be 
characterized by the death of some individuals, with the consequent disappear- 
ance of their surnames, and no significative appearance of new surnames. Since 
the total population grows, however, N{s) will eventually attain a minimum and, 
from then on, will increase. 

To illustrate this situation, we numerically solve eqns (ll) to (ID for the 
initial condition (1,20), with 20 families of one individual each. Note that this 
is the initial condition for which we detected the largest deviations from the 
first-order approximation in Section iA. Figure 5 shows the evolution of N(s) 
for different values of /x. As expected, an initial transient where the surname 
diversity drops is found for /j, > 0. The transient is longer and the minimum in 
N{s) is deeper as fj, grows. In all cases, however, the subsequent growth of N{s) 
is clearly seen. Note that the slope of this growth depends slightly on fi, a feature 
not predicted by the first-order continuous approximation. 

These long transients where (for large fi and suitable initial conditions) the 
number of different surnames is expected to decrease could be relevant to the de- 
scription of modern populations with declining diversity (Cavalli-Sforza & Feld- 
man, 1981; Legay & Vernay, 1999). A realistic description of this situation should 
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Figure 5: Evolution of the average number of different surnames for the initial 
condition (1, 20), with a = 10~^ and different values of the death probability fi. 

however take into account that, in recent times, the relative values of birth and 
death rates in actual populations have considerably changed. The general trend 
to population growth observed worldwide in the nineteenth century has by now 
been reversed in many developed areas, such as in Europe, where the total pop- 
ulation is practically stationary. In our model, this corresponds to an increment 
in the value of to values close to unity. In particular, the death probability 
must be allowed to vary with time. Additionally, if the description is expected 
to encompass the periods where the appearance of new surnames was frequent 
-specifically, the Middle Ages in the case of Europe- a should also change during 
the evolution. 

As a qualitative demonstration of the effect of varying the model parameters 
with time, we focus the attention on a change of fi from a low value to a high 
value. For the initial condition (1,20) we fix /i = during the first sq evolution 
steps. Then, n is left to increase linearly with s, such that it reaches unity after 
si additional steps. In the numerical calculations shown in Fig. 6 we fix si = 10^ 
and consider several values of sq. For s < sq, the number of different surnames 
increases (cf. Fig. 5). Then, as the death probability grows, N{s) attains a 
maximum and begins to decrease. For asymptotically large times, it is expected 
to approach a stationary value, as predicted for the case ^ = 1 in Section |3.3|. 
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Figure 6: Evolution of the diversity of surnames under variation of the death 
probability fi, for a = 10~^. See main text for further explanations. 

Interestingly, N(s) responds faster to the effects of growing /j, for smaller values 
of So, where the contribution of the initial condition is still important during the 
variation of the death probability. The asymptotic surname distribution seems 
to be, in this sense, quite robust to the action of varying n. This feature should 
be related to the fact that the asymptotic distribution is rather insensitive to the 



value of /X, as discussed at the end of Section |3.2. 



5 Comparison of the model with field data 

Finally, we compare some of the analytical results for our model, derived under 



several approximations in Section 3.2, with actual data from three modern pop- 
ulations. Specifically, we focus the attention on the second-order approximation, 
eqn (|2^) , for the asymptotic distribution of surnames in the range of small family 
sizes, where the effect of the initial condition is negligible. In fact, it is virtually 
impossible to extract information on the distribution of surnames in historical 
times from the data presently available. We recall that eqn (^) correctly de- 
scribes the asymptotic power law n{y, z) oc with ~ 2, as observed in real 
data. The validity of our second-order approximation is therefore to be especially 
evaluated in the region of very small family sizes, where the distribution differs 
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from the power law. 

Our three data sets were obtained from surname counts in telephone books. 
They correspond to (i) the almost 350,000 different surnames of the whole 1996 
Argentine telephone book, (ii) the 6,400 surnames beginning by A in the 1996 
Berlin telephone book, and (iii) the surnames in five Japanese cities, with popu- 
lations ranging between 2 x 10^ and 2 x 10^, reported by Miyazima et al. (2000). 
Let us point out that the three populations involved here have considerably differ- 
ent demographic history. The modern Argentine population has predominantly 
European ancestors, who immigrated mainly in the period 1880-1915 and just 
after the World War II. Their surname distribution has therefore to be consid- 
ered as a combined sample from the countries that contributed the immigrants. 
From the times of the largest immigration waves to the present, both the birth 
and death rates have substantially changed. As for Berlin, this city was practi- 
cally abandoned in the late stages of World War II and subsequently repopulated 
with a mixture of the ancient inhabitants and newcomers from other cities of 
Germany. In this case, the surname distribution has consequently resulted from 
a combination of several German regions. The modern population in Berlin, in 
addition, presents the particularity of having been artificially separated into two 
practically immiscible groups during several decades, ending in 1989. Since Euro- 
pean surnames appeared mostly during the Middle Ages (Legay & Vernay, 1999), 
the populations that contributed surnames to both Argentina and Berlin are ex- 
pected to have developed the asymptotic surname distribution on a substantial 
range of family sizes. This situation contrasts with the case of modern Japanese 
surnames, that were originated some 120 years ago (Miyazima et al, 2000). In 
Japan, moreover, the contribution of immigration in the relevant period should 
have been considerably less important than in Berlin and Argentina. 

Explicit evaluation of eqn ( p3|) requires fixing the values of the mutation 
rate a and the death probability fi. We know that the mutation rate, i.e. the 
probability that a new individual acquires a surname not previously present in 
the population, is very low for any modern society. Since, as far as a ^ 0, the 
limit of eqn (^3|) for q ^ is well defined, fixing any sufficiently small value for 
a gives a correct description of n{y,z). In our comparison with real data, we 
have taken a = 10"^. As for the death probability fx, unfortunately, we have 
found it impossible to fix a reliable value. For all real populations, mortality 
has considerably changed during the periods relevant to the evolution of surname 
frequencies. Moreover, since fi measures the relative frequency of death and birth 
events -the latter also including, in real populations, arrival of new individuals by 
immigration- an evaluation of the death probability should also involve a detailed 
description of immigration effects. On the other hand, as discussed at the end 



of Section |3.2| , the asymptotic profile of n(y, z), given by eqn (|23D, is practically 
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Figure 7: Frequency of family names as a function of the family size for three 
modern populations. Dots correspond to real data from the sources quoted in the 
text. Curves stand for the fitting of each data set with eqn (p3|), taking fi = 0.7 
for Argentina and /i = 0.2 for Berlin and Japan. For clarity, the data for Berlin 
and Japan have been shifted by a constant factor. 
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independent of the death probabihty as far as fi is not close to unity. We therefore 
decided to fit the field data with eqn ( p3| ) allowing fj, to vary in order to get the 
best approximation in the relevant zone of small family sizes. 

Figure 7 shows the three data sets and the corresponding fittings with eqn 
(^). Fittings have been optimized in the domain of small family sizes, where 
our analytical approximation is expected to hold. In the case of Argentina, the 
agreement with real data is excellent up to family sizes above i = 100. Only for 
i > 200 a systematic deviation is observed, where the analytical result overesti- 
mates the frequency. According to our discussion in Section |3.2| , this deviation 
would be the remnant contribution of initial conditions. For the Berlin data, the 
statistics are poorer. However, it is clear that the analytical approximation fits 
the data with good precision in the whole range shown here. In the case of the 
Japanese data the fitting is very good for relatively small family sizes, i < 10, 
but, on the other hand, noticeable systematic deviations appear for i > 20. This 
agrees with our expectation on the effect of initial conditions. In the 120 years 
elapsed from the appearance of modern Japanese family names, the population 
in the cities from which the present data were obtained has increased by a factor 
of, at most, order 10. Therefore, according to eqn (p^), initial conditions should 
contribute to the distribution for relatively small family sizes, i > 10, just as 
observed. 

6 Discussion 

We have analysed a birth-death model with overlapping generations, in order to 
study several statistical properties of monoparental inheritance in large popula- 
tions. We have focused the analysis on a cultural trait, namely, the distribution 
of surnames, taking advantage of the availability of big corpora of real data. 
However, the model applies generally to biological traits associated with non re- 
combining alleles. Our model has two parameters, namely, the probability that 
a new individual carries a surname not present in the population, a, and the 
average number of death events per birth /i. For any value of a > the system 
attains a broad, stationary distribution of surname diversity. If /i < 1 the total 
population grows exponentially in time and so does the total number of different 
surnames. The marginal case /i = 1 corresponds to constant total population 
and, on average, to constant surname diversity for asymptotically large times. 

Our analytical results for the stationary distribution of surnames frequency 
are in good agreement with field data for modern human populations in different 
countries. Through an analysis of the transient time required for this distribution 
to reach its asymptotic shape, we have shown that some deviations observed in 
real data might actually reflect the composition of the founder population. This 
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result has implications in the study of polyphyletism. Indeed, if the same sur- 
name can have multiple origins and thus the individuals carrying it are not always 

phylogenetically related, this will affect the shape of the surname distribution. 
In particular, it is not difficult to estimate the time when surnames originated 
in a population (using historical records) or, in the biological counterpart, when 
a mutant allele first appeared (through molecular clock analysis). Then, the ap- 
proximate cut-off until which the stationary distribution follows the asymptotic 
shape can be calculated. If the distribution underestimates the frequency of val- 
ues larger than the cut-off, the system is mainly polyphylctic. If it overestimates 
that frequency, then the simultaneous appearance of many individual carrying 
different surnames took place in the past. 

The strong resemblance between the cultural inheritance of the surname and 
the biological process in which non recombining neutral alleles are passed to 
offspring has justified to apply results from field data in the former case to the 
latter (Barrai et al., 1996). In the few cases where data on genetic diversity was 
available, it was possible to retrieve information on past populations by comparing 
both sets of data (Sykes & Irven, 2000). A specific example comes from the small 
island of Tristan da Cunha, where the fact that its 300 inhabitants represent only 
seven surnames and five mitochondrial lineages reflects without doubt the small 
size of the founder population (Soodyall et ai, 1997). 

Our results could be applied to population genetics under some hypotheses. 
Indeed, we are assuming that the number of different alleles at a given locus 
is practically infinite (this is analogous to the assumption made by Kimura & 
Crow (1964), in their model of infinitely many alleles), since the possibility of 
backward mutations is discarded. Nonetheless, this factor could be accounted 
for just by lowering the value of a, because this process implies that the "new" 
mutant is in fact identical to one of the forms present in the population. Since, 
as long as a is small, our results are not sensibly modified, we could also work 
with a finite but large number R of different surnames (equivalently, large genetic 
polymorphism), and use the same model as long as the current diversity is lower 
than R. We have also assumed that the death rate is constant during the life time 
of individuals, while it is known that the life expectancy depends not only on the 
age of individuals, but also varies as a function of time (Vaupel et al, 1998). A 
more realistic model of human inheritance could be constructed by taking into 
account the variation of // along the lifetime of each individual. 

The evolution of surname diversity follows truly neutral evolution: family 
names do not fulfill any practical purpose but identifying the lineage of each 
individual. They cannot be selected for through natural selection. Hence, as 
we have shown, their statistical distribution closely follows the predictions of a 
neutral model of monoparental inheritance. It would be interesting to test our 
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theoretical results against the genetic diversity of a large population sample. Un- 
fortunately, however, data on frequencies of allelic polymorphisms are still scarce 

to carry out a massive study like the one presented here for actual surname dis- 
tributions. Moreover, one should make sure that the different haplotypcs in the 
sample do not confer any advantage to the individuals, in which case positive feed- 
backs and deviations from neutral statistics would be expected. Different tests 
have been proposed in the literature to detect deviations from neutrality (Nielsen, 
2001; Fu & Li, 1993), and a number of neutral haplotypes have been positively 
identified (Sanchez Mazas et al, 1994; Stcnoicn, 1999). Hopefully, enough data 
will be available in a near future to calculate reliable diversity distributions in 
human populations. 
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A Average evolution of the total population 

As explained in the main text, at each step, the total population P{s) in our 
model either increases by 1 with probability 1 — /U, due to the occurrence of a 
birth event and no death event, or remains constant with probability due to 
the occurrence of both a birth and a death event. Equations (|^) and (^) quantify 
the corresponding stochastic process. 

Since the total population changes during the evolution, the time interval 
6t to be associated with each step -which corresponds to the interval between 
consecutive births- must also change. In fact, the birth frequency is proportional 
to the total population, so that 5t is inversely proportional to P{t). We write, at 
step s. 



where the frequency fixes time units. 

In consequence, the real-time average variation in P{s) can be obtained from 



where overlines indicate average over realizations of the stochastic process w{s). 
To obtain this average evolution equation we have used eqn (||), and have taken 
into account that w{s) and P{s) are independent stochastic processes at each 
step s, so that w{s)P{s) =w{s)P{s) = fiP{s). 

If v is identified with the birth rate per individual and unit time, the product 
ufj, is the mortality rate. For constant and fi, the solution to eqn ( ^2| ) is 




(31) 



dP _ lP{s + l) - P{sy 
Ht ^ [ 6t{.s) . 



u[l - w{s)]P{s) = z/(l - n)P, 



(32) 



P(t) = Poexp[z/(l-/i)t] 



(33) 



with Pq the initial population. 



27 



B Continuous approximation to the distribution of 
family sizes 

Under the action of both birth and death events, the evohition of the average 
number of famihes of size i, Hi, is given by eqns (|l3|) and (^). These equations 
can be approximately solved assuming that the solution varies slowly on s and i, so 
that these two discrete variables admit replacement by continuous variables z and 
y, respectively. Accordingly, ni{s) is replaced by a continuous function n{y,z). 
The approximation is based on the expansions ni±i{s) = n{y, z) ± dyn{y, z) + • • • 
and ni{s + 1) = n{y, z) + dzn{y, z) + • • •, to be introduced in eqns ([l^ ) and (|l5|) 
at different truncation orders. 



B.l First-order approximation 

To the first order, we obtain the differential equation 

dz _Po + (1 - oy 

for y > 1. The contribution for y = 1, given by eqn ([l^, can be incorporated 
to (|3^) either as a boundary condition or as a singular inhomogeneity in the 
equation itself. The latter yields 

c)tI 1 — Q! — d 

li~ + ^~TTi ^T7-{yn) = a5{y - I), (35) 

oz Po + (1 - oy 

where 5{x) stands for the Dirac delta distribution. In the following, we obtain 
and analyse the solution to this equation. 
Introducing the auxiliary variables 

C = ln('l + i^zV ?7 = lny, (36) 



eqn (|35D becomes 



a4 1 — jj. or] I — fi 

with f{i],(,) = yn(y,z). This is a one-dimensional wave equation, with "spatial" 
variable rj and "temporal" variable It describes shape-preserving advection of 
the "density" / with constant velocity v = {1 — a — — fx), subject to the 

action of a point source of intensity a/(l — /i) at r/ = 0. In our problem, thus, 
this equation makes sense for f > (/i < 1 — a) since rj must be non-negative 
(y = ^ ^ 1)- The general solution for 77 7^ is given by an arbitrary combination 
of functions of the form f{i],(,) = fo{i] — v^), where /o is in principle arbitrary. 
The combination must be chosen in such a way that the boundary and initial 
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conditions are satisfied. In our case, this is achieved with a combination of two 
such functions. In the original variables, the solution is given piecewise as 

1 — a — fi 

for y < yriz), and 

n{y, z) = y:^^n{y/yT{z), 0) (39) 

for y > yriz). Here, n(y, 0) is the initial distribution. The transition point 
between the two pieces is located at 

yT{z)=(l + —^z) . (40) 



These expressions give the first-order continuous approximation n{y, t) to the 
distribution of families by size. 

The average total number of surnames in the continuous approximation is 
defined in (21). For the first-order approximation we find 

N{z) = r n{y, z)dy = + (1 - /^)^ ^-c^^ + T ^(^/^^^ q)^. (41) 
Ji 1 - a — fj, Ji Jyj. yT 

The last integral is clearly equal to the initial number of surnames, Nq, as may be 
immediately realized by the change of variables y/yx U- Explicit calculation 
of the first integral shows that N{z) = Nq + az, i.e. exactly the same result as 
for the case of /i = 0, eqn ([Tl|). As a function of time, we have 



N{t) = No + -^{exp[z.(l - fj,)t] - 1}. (42) 



B.2 Second-order approximation 



Truncation to the second order in the continuous approximation to eqns (|13| ) and 
(0) yields 

dz Po + (1 - l^)z dy 2[Po + (1 - fi)z\ dy^ 

Unfortunately, this equation cannot be analytically solved for arbitrary initial 
conditions. However, a particular solution can be found in the form of a separate 
function, n{y,z) = h{y)P{z) = /i(y)[Po + (1 ~ t^z]. Comparison with eqns @ 
and (1£) suggests that this particular solution will correspond to the long-time 
asymptotic evolution. It reads 

nfe..) = -r^^^ U}-^f\-^uU- 1,0, 2l^,) . (44) 

l — a — ^\ 1 — a + fiJ \ 1 — a + fi J 

where U (a, b, x) is the logarithmic Kummer's function (Abramowitz & Stegun, 
1970). 
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